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Abstract 

The neutrino energy density of the Universe can be conveniently parametrized in 
terms of the so-called effective number of neutrinos, N^^ . This parameter enters 
in several cosmological observables. In particular it is an important input in those 
numerical codes, like CMBFAST, which are used to study the Cosmic Microwave 
Background anisotropy spectrum. By studying the neutrino decoupling with Boltz- 
mann equations, one can show that this quantity differs from the number of massless 
neutrino species for an additional contribution due to a partial heating of neutrinos 
during the annihilations, leading to non thermal features in their final distribu- 
tions. In this paper we review the different results obtained in the literature and 
perform a new analysis which takes into account, in a fully consistent way, the QED 
corrections at finite temperature to the photon and plasma equation of state. 
The value found for three massless active neutrinos is N^^ = 3.0395, in perfect 
agreement with the recommended value used in CMBFAST, N^^ = 3.04. We also 
discuss the case of additional relativistic relics and massive active neutrinos. 
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1 Introduction 



At temperatures below the muon mass and above ~ 10 MeV, the Universe is 
filled by a plasma of photons, electrons, positrons, and neutrinos, kept in ther- 
modynamical equilibrium by the electroweak interactions. As the temperature 
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drops below this value, the rate of weak interactions starts to be comparable 
with the Universe expansion rate, and neutrinos decouple from the electro- 
magnetic 7, plasma. For most practical purposes, it is accurate enough to 
consider the freeze-out of neutrinos as fully achieved at temperatures of about 
2 — 3 MeV. In this limit neutrinos do not share any entropy release from 
annihilations, once the temperature drops further, below the electron mass. 
Assuming that all the entropy produced by the annihilations is transferred to 
photons, their temperature T is increased with respect to the neutrino tem- 
perature by the well known factor T/T^ = (11/4)^^'^ (see, for example, [1]). 
Actually more accurate calculations show that neutrinos are still slightly in- 
teracting with the electromagnetic plasma and thus receive a small portion 
of the entropy from e"^ annihilations [2]. Neutrinos with higher momenta are 
more heated, since, in the relevant range of energies, weak interactions get 
stronger with rising energy. This produces a momentum dependent distortion 
in the neutrino spectra from the equilibrium Fermi-Dirac behaviour. 

A further, though smaller, effect on T /Ty is induced by finite temperature QED 
corrections to the electromagnetic plasma. In fact electromagnetic interactions 
modify the and 7 dispersion relations, and thus the energy density and 
pressure of the plasma. More precisely, the energy density is lowered so that the 
annihilation phase releases less entropy with respect to the non-interacting 
particle limit calculation. Since most of this energy ends up into photons, this 
decrease results in a smaller T/Ty ratio [3]. 

While any direct observation of the actual distortions in the neutrino distri- 
butions is out of question, nevertheless we can hope that the increase of the 
total energy of the relic neutrinos may have a sizeable effect on the expansion 
rate of the Universe. The energy density stored in relativistic species, pR, is 
customarily given in terms of the so-called effective number of neutrino species 
[4,5], Nl^, through the relation 



where is the energy density of photons, whose value today is known from 
the measurement of the Cosmic Microwave Background (CMB) temperature. 
Eq. (1) can be also written as 



where denotes the energy density of a single specie of massless neutrino with 
an equilibrium Fermi-Dirac distribution, and p° is the photon energy density 
in the instantaneous neutrino decoupling approximation. The normalization 
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of A^^ is such that it gives A^^ = 3 in the standard case of three famihes of 
massless neutrinos, in the hmit of an instantaneous decouphng. In principle, 
A*"^^ can receive contribution from other relativistic rehcs with energy density 
Px as well. In the following we will mostly restrict our analysis to the standard 
case, but we will further consider this more general framework in Section 4. 

As we will discuss below, when considered separately, the non instantaneous 
neutrino decoupling gives AA"^^ = A^^''^ — 3 ~ 0.034 [6-8], while QED effects 
contribute for about AA^^^ ~ 0.011 [3]. Could the two effects be added linearly, 
they therefore would produce a final value A"^^ ~ 3.045. Recently, these effects 
have been reconsidered in [9]. This work combines the results for the non 
instantaneous neutrino decoupling obtained by a numerical calculation in [10] 
and then replacing the ratio T/Tj, = (11/4)^''^ with the value obtained by 
considering QED corrections. This procedure may provide a good estimate, 
but it is worth pointing out that, since QED effects modify several points of the 
Boltzmann equations describing neutrino decoupling, it is advisable to study 
any possible interplay of the two effects. Actually a precision calculation of this 
kind, where both these effects are included at the same time, is still lacking, 
and it is the aim of the present paper. We will show that this interplay ends 
up in a 10% smaller total correction to A"^^ than what would be obtained by 
simply adding the two contributions. We also notice that the calculation of [10] 
seems to underestimate^ the corrections from non instantaneous neutrino 
decoupling with respect to the result found in [6-8]. Our analysis is based 
on a numerical code described in [8], which has been modified to take into 
account QED finite temperature effects. 

From the observational point of view the effects considered in this paper are 
too small to influence Big Bang Nucleosynthesis (BBN), since they produce a 
change in the ^He mass fraction, Ay(^He) ~ 10~^ [8,11,12], which is smaller 
than the actual theoretical, 5 10~^ [12,13], and experimental, > 2 10^'^ [14], 
uncertainties on this quantity. This change is actually slightly larger when one 
takes into account the effect of flavour neutrino oscillations, as recently shown 
in [15]. More promisingly, they might be detected via future precision CMB 
anisotropy measurements at high multipoles, since cosmic variance prevents 
their resolution on scales probed by present balloon experiments. According to 
a recent analysis [16] (see also [13,17-21]), CMB temperature measurements 
by the Planck satellite experiment will be able to provide a measure of the 
relativistic energy density with a precision of about AA^;;^ ~ 0.2 , without 
any strong prior on the other cosmic parameters. Furthermore, as discussed 
in [22], the situation is more promising if polarization measurements will be 



^ In ref. [10] it is claimed that their result is enhanced with respect to the one 
of [6] . We think that this is actually due to a misinterpretation of the results of [6] 
in terms of A'^^^ . This can be easily understood by noticing that in [10] it is found 
a final ratio T/T^ closer to (11/4)^''^ than what obtained in [6-8]. 
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available, and some stronger priors are imposed on other cosmological param- 
eters, enforced for example by independent measurements. In this case, AA''^''' 
would be determined with an accuracy comparable or, according to [22], even 
higher than the order of magnitude of the effects we are here considering. We 
notice that in the CMB data analysis the presence of a non vanishing AN^^ is 
already considered with, for example, a recommended default value for three 
massless active neutrinos A^^''^ = 3.04 in the CMBFAST code [23] . A precise 
check of this important parameter is therefore mandatory. 

The paper is organized as follows. In Section 2 we summarize the set of equa- 
tions and our numerical approach to compute the neutrino distortion due to 
their incomplete decoupling during the annihilation phase. In Section 3 
we also consider the corrections to the equation of state of and 7 plasma 
due to QED effects at finite temperature, showing how these effects modify 
the numerical computation of Section 2, and we present our results. We also 
discuss the general case of extra relativistic degrees of freedom and of massive 
active neutrinos in Section 4. Finally our conclusions are reported in Section 



2 Corrections from to the non instantaneous decoupling 

The effects of non-instantaneous neutrino decoupling have been addressed in 
several studies, either based on analytical methods [2], on numerical analy- 
sis with some simplifying approximations [24-27], or finally on full numerical 
computations, which solve the Boltzmann equations until neutrinos are com- 
pletely decoupled [6-8,10,28]. 

The procedure is easily summarized. Following the notation of [6,8], a con- 
venient choice for the time variable is given hj x = niR, where R{t) is the 
scale factor of the Universe (chosen to have dimension of length) and m some 
reference mass, taken to be the electron mass as in [8]. We also introduce the 
dimensionless comoving momentum y = pR, and the rescaled photon tem- 
perature z = T R. In this notation, the Boltzmann equations for the neutrino 
distributions can be written in the form 



where H is the Hubble parameter, while 1,^^ represents the coUisional integral, 
in momentum space y, for the single neutrino specie z/^, and is a functional of 
all neutrinos and electron/positron distributions^. For the processes of interest 

^ Since are kept in thermodynamical equilibrium with 7 by the electromagnetic 
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(see [6] for the detailed calculation), some of the integrations appearing in the 
Iij^ can be analytically performed, and the collisional integrals can be reduced 
to two-dimensional integrals. In the range of temperatures we are interested in, 
electron neutrinos experience charged current interactions due to the presence 
of in the thermal bath, while all active neutrino flavours interact via neutral 
current interactions. For this reason, the noncquilibrium corrections to the 
distribution function fy^ are different from those of the other two neutrino 
species fy^ {x denoting both fi and r). 

The two Boltzmann equations for f^,^ and f^^ must be supplemented by the 
continuity equation 

A,W = i(,-3P), (4) 



where p = p{x/me) , P — P{x/me) , and p and P are the energy density 
and pressure of the 7, e^, u plasma. Eq. (4) can be rewritten as an evolution 

equation for the quantity z which gives the ratio T/Ty. In case of instantaneous 
neutrino decoupling, the asymptotic value of z, denoted by z^ , results to be 
the well-known value (11/4)^/^. 

As in Ref. [8] the unknown neutrino distributions are parametrized as 



1 



1 + at{x) P, {y) 

i=0 



(5) 



where Pi{y) are orthonormal polynomials with respect to the Fermi function 
weight 



By substituting Eq. (5) into Eqs. (3), one gets 
, 1 °° 

Tx ''^^''^ I'^y' ^y'^ f^'^- ■ 


Since at high temperature the neutrinos are in thermal equilibrium, the initial 
condition for the coefficients is of = 0. 

interactions, they have Fermi-Dirac distributions and share the same temperature 
of photons, T. We neglect the completely irrelevant asymmetry. 
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Flavour {a) 


Co 


cf 




^3 


e 


-2.556 


-2.739 


6.133 


-0.1477 


11, T 


-2.049 


-2.259 


3.145 


-0.1129 



Table 1 

Values of the coefficients of Eq. (9) from [8] (QED effects not included). 

To perform the numerical computation it is necessary to truncate the infinite 
series in Eq. (5) at some finite value n, where the choice of n depends on the 
accuracy we do require on the results; it is found in [8] that n = 3 gives an 
accuracy of about 1% in the neutrino distortions. In this case, the asymp- 
totic value found for z, denoted by z^"^, slightly differs from the instantaneous 
decoupling result z^ 

^ ^ = -1.406.10- . (8) 

The final neutrino distribution functions show a non thermal behaviour due 
to the presence of non vanishing ripple terms. Having fixed n = 3, we rewrite 
the expression of Eq. (5) in the form 

and report the final values for the coefficients in Table 1. By using Eq. (9) 
it is now immediate to compute the additional contribution to the neutrino 
energy density due to incomplete decoupling. Defining 



'° J,°° dy y^ (ev + ly' 



= 10- (^c^ + 2700 [y^jct + 310 [Y^jc^2+ 12150 (^^J ^3 J , 

(10) 



one obtains, for the values of Table 1, Sp^^/ p° = 0.953% and Sp^y^/p^ = O.oyy 7o. 
Finally, from the definition of Nf^ of Eq. (1), it is straightforward to get the 
following expression 



.^'"7 VP. Pu J \ p- p\ 



V 



which gives the effective number for three massless active neutrinos for tem- 
peratures below the neutrino decoupling (~ 1 MeV). The values found from 
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the mimcrical analysis then give = 3.0345, which fixes the additional con- 
tribution to be AiV;;'^ = 0.0345. Notice that from Eq. (11), ANf takes two 
contributions. The term proportional to 6z/z^ accounts for the smaller profit 
that photon temperature gets from e"^ annihilations, whereas the contribu- 
tions due to Spu^ ^/pl are a measure of the non thermal behaviour of neutrino 
distribution functions. These two terms are indeed of the same size. 

Remarkably, our results are in very good agreement with a previous analysis 
performed applying a different numerical technique [6,7]. They find Sz/z'^ — 
(— 1.37±0.02)-10~^, and an increase of energy with respect to the instanta- 
neous decoupling case of Sp^^/pl = (0.946±0.001) % for the electron neutrinos 
and of Sp^Jpl = (0.398±0.001) % for the other two species. These values give 
the effective number of neutrino species Nf^ = 3.0340±0.0003 . 

On the other hand we disagree with the results of [10], where by solving the 
same Boltzmann equations it is found Sp^^/pl = 0.607%, Sp^^/pl = 0.256% 
and 5z/z° = -0.888-10-3, which finally gives Nf = 3.022. In both calcula- 
tions [6,7] and [10], the integral-differential Boltzmann equations are solved 
through a grid in momentum space. The grid adopted in [10] extends to a 
larger momentum range than the one of [6,7]. The grid used in [6,7] is how- 
ever denser in the interval 0.1 < y (m/MeV) < 20 . For neutrino distributions 
in thermal equilibrium, the 97.5% of the energy density comes from particles 
with momentum in the range 1 < y (m/MeV) < 10 . According to the analy- 
sis of [7], the different results of [10] may be due to a lack of precision in this 
most relevant interval. To conclude, we can only point out that the beautiful 
agreement of our findings with what obtained in [6,7], despite of the rather 
different numerical method, make us rather confident on our result for Nf^. 



3 Corrections from QED at finite temperature 

Finite temperature QED corrections modify in several points the calculations 
of neutrino decoupling described till now. First, through the change on the 
electromagnetic plasma equation of state, they affect Eq. (4). Moreover, since 
masses are renormalized by a finite temperature term, the r.h.s. of Eq.s (3) 
and (7) should be modified correspondingly. Finally, the change of energy den- 
sity modifies the expansion rate H in the Boltzmann equations (3). The effects 
of these corrections on neutrino decoupling temperature have been considered 
in [29]. Their analysis is however in the framework of the instantaneous neu- 
trino decoupling limit, therefore they do not consider the non thermal effects 
described in the previous Section. 

The change in the electromagnetic plasma equation of state can be evaluated 
by first considering the corrections induced on the and photon masses. 
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They can be obtained perturbatively by computing the loop corrections to 
the self-energy of these particles. For the electron/positron mass, up to order 
a = e^/ (47r) we find the additional finite temperature contribution [3] 



O 71 



dk 



e 1 

E~k e^fc/^ + 1 



2 ml ^ 
up 



dk — log 



p + k 



p 



k 



1 



'.m 



where = k"^ + ml- While the first two terms of this expression depend on 
the plasma temperature T only, the last one depends on the momentum p 
as well. However, by averaging this term over the equilibrium distribution, 
one easily finds that it contributes for less than 10 % to 5ml ■ -^^^ ^^^^ reason 
we will neglect it in the following (see also [12]). 

The renormalized photon mass in the electromagnetic plasma is instead given, 
up to order a, by [29] 



5m' (T) 



dk 



Ek e^fe/^ + 1 ■ 



(13) 



The corrections (12) and (13) modify the corresponding dispersion relations 
as Ef = k"^ + mf + 5ml (T) {i = e, 7). This in turn affects the total pressure 
and the energy density of the electromagnetic plasma 



m 

p= fdkk' log 
J 



P- 





dP 

-P + T — 
dT 



(1- 



-Ee/ry 

-E-f/T^ 



(14) 



(15) 



Expanding P with respect to 5ml 5m^, one obtains the first order cor- 
rection 



P 



int 



dk 

2^ 



k 5m'^ (T) 
Ek e^^l'^ + 1^2 e^l^ - 1 



5ml{T) 



(16) 



Note that an important factor 1/2 must be introduced for not double counting 
the renormalization effect on the total pressure, which now reads P = po+p™* 



^ Our result agrees with Eq. (B2) reported in [3], while we think that the analogous 
result in Eq. (35) of [12] contains some misprints. 
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Flavour {a) 


Co 


cf 




^3 


e 


-2.507 


-2.731 


6.010 


-0.1419 


11, T 


-2.003 


-2.196 


3.061 


-0.1091 



Table 2 

Values of the coefficients of Eq. (9) when QED effects are included. 

[3], being the value of the pressure for the non- interacting particle gas. The 
energy density is then obtained by using P in Eq. (15). 

The presence of the additional contributions P'"' and p'"*^ modify the evolution 
equation for z contained in Eq. (4), which now reads 

^ ^ ^A^h) - 2^^ r dyy^ + 2^) + G,{x/z) 

where 



Gi{uj) = 2na 



\ D D 



(18) 

+2™ a; 1^^^^ - K{uj) K'{uj) + + J'{uj) K{uj) + J{uj) K'{uj)^ , 

(19) 



with 



oo 2 



1 /" U 

K(uj) = — du , ^ , (20) 

{ exp (Vi?+^) + 1 ' ^ ^ 



du ^ A — , (211 



1 ^ 

-I "" " (exp + cj2) + l) 



1 y ^ exp ( + a;2 ) 
^ i (exp + l) 



The functions K'^uj) and J'{u!) stand for the first derivative of K^uj) and J(a;) 
with respect to their argument. Note that neutrinos affect the final value of z 
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[D, < J (^no-Vr^ii/jj ) 


[iuj i^no-v^ii/JJ ) 




our worK (^Vc^H/iJ j 


ozj z 


( 1 o7_Ln ^o^1^— 3 
\ — \..6l±\i.\il)\\} 


— U.ooo-lU 


1 /in/^ 1 n— 3 
— 1.4U0-1U 


— i.o4i- iU 




(0.946±0.001)% 


0.607% 


0.953% 


0.935% 


^PuJpI 


(0.398±0.001)% 


0.256% 


0.399% 


0.390% 


^ ^ V 


3.0340±0.0003 


3.022 


3.0345 


3.0395 



Table 3 



The results of the different analyses. 

through the terms dfu^/dx, which are not vanishing only if neutrinos have a 
non thermal behaviour. 

In case one neglects the finite temperature QED corrections, the functions 
Gi{x/ z) and G2{x/ z) in (17) vanish, and one recovers the expression reported 
in [8]. Notice that the presence of G2{x/ z) in the denominator of the r.h.s. 
of (17) makes, at least in principle, not correct to simply sum the neutrino 
contribution with the QED one. 

Once both effects, neutrino incomplete decoupling and QED corrections to 
electromagnetic plasma, are included into the code of Ref.[8], we find the new 
results Sp^Jpl = 0.935%, Sp^Jpl = 0.390% and 5z/z^ = -1.841-10-3 which 
give Nl^ — 3.0395. The values of the cf coefficients can be found in Table 2. In 
Table 3 we report a comprehensive summary of our results. Comparing the last 
two columns of this table, we see that introducing the finite temperature QED 
corrections in the non instantaneous decoupling scenario, leads to a change 
on the effective number of neutrinos of 0.005 , which is a factor 2 less than 
what has been obtained in Refs. [3,12]. This is actually due to the interplay 
of incomplete decoupling of neutrinos and plasma effect (see the r.h.s of Eq. 
(17)), which therefore, at the level of accuracy we are considering, cannot 
be considered separately and then added linearly. Notice that the increase 
of Nf^ is mainly due to a variation of 5z/z^, the changes on 5py^/ p^ being 
much smaller. Remarkably, the overall result Nl^ — 3.0395 turns out to be in 
excellent agreement with the recommended default value Nl^ — 3.04 used in 
the CMBFAST code [23]. 



4 Extra relativistic degrees of freedom and massive neutrinos 

We now consider the possibility that, at the stage of neutrino decoupling, there 
are extra relativistic degrees of freedom, provided by some X field excitations, 
which are assumed to have a thermal distribution with some temperature 
Tx- Their contribution px to the total energy density of the Universe can 
be parametrized in terms of an additional contribution AA^^^ in the effective 
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number of neutrinos, as defined in eq. (2) 

={^] Nx, (23) 




where 



11x4/3 



and we have defined zx — TxR- The parameter gx depends on the spin 

{9x = 1 for a real scalar, gx = 7/4 for a Weyl spinor, etc.) as well as on the 
additional internal degrees of freedom of the X particles. Notice that if the X 
excitations have decoupled between /i"^ and annihilation phases we simply 
have Nx — 4/1 gx- For an earher decoupling we have instead Nx < 4/7 gx- 

The presence of px, apart from introducing a new contribution to Nf^ , slightly 
affects the results obtained in the previous Sections, namely the relative change 
of neutrino energy density 5py^ ^/p^ induced by incomplete neutrino decou- 
pling, as well as the asymptotic photon temperature z^"^. In fact since their 
energy density increases the expansion rate of the Universe, we expect 5p^^ ^/ p° 
to decrease with growing Nx-, and the ratio z^ / z^^ to become closer to unity, 
since neutrino decoupling process starts at earlier times. If we denote with 
5pue,SNx) / p^ and z^^{Nx) the new values for these parameters as functions 
of Nx we therefore have 



Since we are interested to those contributions to Nl^ corresponding to species 
which are relativistic for temperatures in the MeV range, we can severely 
bound /\Nf using the results from BBN, which leads to the conservative 
bound /\Nf < 1 ( see for example [30] for a recent analysis). Using our nu- 
merical code we have evaluated how z^^{Nx)/z^ and Sp^^{Nx)/p^ change 
when Nx varies in the range < Nx < 1. In this interval, with an accuracy 
of 10-^ on Nf, we find 



3.0395 + (1 - /3) 




+ (1 - P)Nx 



(26) 
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with f3 ~ 0.0014 and where now both j z^^ and ^Pv^^l P% this expression 
are those reported in the last column of Table 3, i.e. with Nx = 0. The 
changes of these parameters with Nx is encoded in an additional term, which 
is weighted by the small parameter /3. 

It should be clear from what we said before that only species which are rel- 
ativistic at the neutrino decoupling, down to the annihilation stage, con- 
tributes for this extra term. It is known that the effective number of neutrinos 
at later stages, as for example at recombination, is much less constrained 
from data [16-21,30,31] and it is still possible the case that energy density in 
the form of relativistic plasma may be injected only well after the Big Bang 
Nucleosynthesis epoch. This implies that the value of A^"^^ may well be rather 
different at the CMB and BBN epochs. In this case, Eq. (26) at recombination 
reads 

Nf = 3.0395 -{1-/3) (^] Ar|s^ + AA^^mb ^ (27) 



where Nx^^ and ANx'^^ are the contribution of species which are relativistic 
at the BBN and neutrino decoupling, and recombination epochs, respectively. 

We finally consider the case of massive active neutrinos. This nowadays plau- 
sible scenario affects, in general, both the expectations for CMB anisotropy 
and large scale structure formation. As it is seems more and more clear from 
neutrino experiments on solar and atmospheric neutrinos, it is unlikely that 
neutrino mass differences may be greater than ~ 0.1 eV [32], unless we enlarge 
the standard scenario introducing sterile neutrino states. At the same time, it 
is also quite well established from Tritium decay data that z/g mass is bound 
to be smaller or, at most, of the order of 1 eV. It is therefore clear that in 
this scenario all neutrino masses are completely negligible as far as their de- 
coupling is concerned. However if their values is as large as ~ 1 eV, they start 
to be relevant as the temperature approaches the range relevant for CMB, 
and the presence of a finite mass modify of course the neutrino contribution 
to A"^^. It is interesting to consider how the effects of incomplete decoupling 
and QED thermal effects studied in the previous Section would now affect 
the neutrino energy density. This is conveniently parametrized by the (time 
dependent) quantity 

Spu^ (m^) ^ p^^jm^) - pl^jma) 

which is defined as in Section 2, but for a neutrino with a finite mass 771^. 
Since for rria ~ 1 eV active neutrino are fully relativistic for temperatures of 
the order of MeV, so we can still take at the annihilation phase the results 
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Fig. 1. Evolution of ^Pu^ 

{nia) / p^y^inia) , for a neutrino mass ma = 1 eV (see text 

for further details). 

of Section 3, it is easy to see that 6p^^{m)/ p1^{m) is given by 



dy ^y^ + (^) (e^ + l)"^ ELo cfy' 
10 , _ = . (29) 



In Figure 1, we plot 6p,y^{ma)/ p^^{ma) for the reference choice nia = 1 eV, 
using the coefficients c" of Table 2. The x range corresponds to a variation 
of the temperature from T ~ 1 keV till T ~ 10^'^ eV. We see that the effect 
of incomplete decoupling and QED plasma masses for and 7 on neutrino 
energy density decreases as the temperature becomes comparable with the 
value of the neutrino mass. Notice that, from (29), at very low values of T, 
when the neutrino energy is dominated by the mass term, Sp^^^nia) / p1^{ma) 
reaches an asymptotic value, which is given by the change in the neutrino 
number density due to the above effects, normalized with the number density 
for a pure thermal Fermi-Dirac distribution. 



5 Conclusion 



In this paper we have considered how the two effects due to incomplete neu- 
trino decoupling and QED corrections to the and photon plasma equation 
of state affect the effective neutrino degrees of freedom, a crucial parameter 
for many cosmological observables. 
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The main result of our analysis, which have been carried out by numerically 
solving the set of Boltzmann equations describing neutrino decoupling, is a 
value for Nf^ = 3.0395, for the standard case of three active neutrino flavours. 
The issue is certainly not new, but we hope our study will contribute to reach a 
better accuracy in the theoretical determination of relativistic energy density 
of the Universe. In particular we have stressed that there is quite an interesting 
interplay between the two effects we have considered, so that at the level of 
accuracy of 0.005 on Nf^, their corrections to the neutrino energy density 
cannot be naively summed as they were fully independent. 

We have also considered the less standard scenario of extra species contributing 
to A^^''^ , and how they affect the calculation of the thermal distortion of 
neutrino distribution, as well as the final value for photon temperature after 
the annihilation phase. 

Despite of the smallness of the corrections we have been concerned with in this 
paper, nevertheless a careful analysis of data on the CMB anisotropy spectrum 
at large multipoles may be able in the near future to reveal their effects. Esti- 
mated sensitivity of Planck satellite experiment to relativistic energy density is 
at least of the order of AA^^^ ~ 0.2, but it is strongly improved, AA^^^ ~ 0.005 
[22], when including polarization measurements and strong priors. 
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